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ABSTRACT 

Finite difference approximation of transonic flow problems is a well- 
developed and largely successful approach. Nevertheless, there is still a 
real need to develop finite element methods for applications arising from 
fluid-structure interactions and problems with complicated boundaries. In 
this paper we introduce a least squares based finite element scheme. It is 
shown that, if suitably formulated, such an approach can lead to physically 
meaningful results. Bottlenecks that arise from such schemes are also 
discussed. 
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1. INTRODUCTION 


In this paper we consider the approximation of transonic flows by finite 
element methods based on a variational method of the least squares type. The 
objective here is purely computational. In particular, we have sought to 
fully exploit the ideas arising from mathematical analysis of such methods 
(see, for example, [1] - [6]) and directly apply them to a nontrivial 
transonic flow problem. The major conclusion drawn from this work is that 
finite element methods — suitably formulated — can give physically meaningful 
results. 

There is a significant and largely successful array of finite difference 
techniques for transonic flows (e.g., [17]). Nevertheless, an assumption 
implicit in this work is that there is still a need for stable and accurate 
finite element approaches. First, there are applications from fluid-structure 
interactions that would benefit from the availability of a finite element flow 
model. Second, there is the issue of complicated boundaries in the flow 
field. The importance of the finite element ideas in such a context — while 
largely untested — is still promising. 

Variational principles of the least squares types have a number of 
valuable computational properties. For example, the algebraic system 
generated is always Hermitian semidef inite. In addition, such schemes, if 
properly formulated, are insensitive to equation type, be it hyperbolic 
(supersonic flows) or elliptic (subsonic flows). In fact, the majority of the 
finite element ideas that have been used for hyperbolic problems to date tend 
to be either implicitly or explicitly of the least squares type. 

Least squares based schemes do have, however, some major computational 
defects. First, they tend to be sensitive to singularities and 
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discontinuities in the flow variables. Moreover, mesh refinement alone does 
not overcome these defects [7]. Based on the work in [7] we introduce 
weighted least squares variational principles, which in combination with mesh 
refinement is capable of dealing with shocks in the flow field. 

In Section 2 we describe the basic numerical formulation, and outline the 
essential computational properties associated with the approach. A key 
feature is the proper choice of weighting functions to use in the least 
squares functional. A closely allied issue is the density modifications 
needed to rule out nonphysical expansion shocks. 

In Section 3 we present sample numerical results. As a model problem we 
select the planar potential flow over a cylinder. 

Other authors have considered finite element approximation of transonic 
flows. Selected references are [18] - [21]. 


2. THE LEAST SQUARES FORMULATION 

A 

We consider the potential flow over a body fi. Let u_ denote the 
velocity and p the density. Then a mass balance yields 


In addition, we have 


div[pu] = 0. 


u = grad 


( 2 . 1 ) 


( 2 . 2 ) 


for the velocity potential 4> • The density p is given as a function of _u 
by the Bernoulli equation. The system is closed by specifying the normal 


velocity 
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_u*n_ = v (2.3) 

A 

at the boundaries of the flow region. On the body 0 the no flow condition 

u_- n. = 0 

applies. We assume that the flow region is contained in a box B and that 
(2.3) is specified on the boundary of B. Thus 

A 

ft = B/n (2.4) 

defines the flow region, and (2.1) - (2.2) hold in ft with (2.3) holding on 
the boundary T and ft. 

Since the flow is assumed to be irrotational , (2.1) - (2.2) can be 

replaced with 

div(pu) = 0 in ft (2.5) 

curl(u) = 0 in ft (2.6) 

u*n = v on T. (2.7) 

A least squares scheme based on this system takes the form 

/ {|div(p_u)|2 + |curl(u)|^} = min, (2.8) 

ft 

where the variation is taken u in some finite element space satisfying the 
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boundary conditions (2.7). Such a div - curl system has proven to be very 
effective for elliptic systems (subsonic flows) in cases where the density 
p = p(u) and the velocity field ii are smooth [8]. 

Preliminary results indicate that with appropriate weighting functions on 
the terms in (2.8), the nonsmooth cases can be treated as well. Nevertheless, 
in this paper we shall focus attention on (2.1) - (2.2) and least squares 
schemes of the form 


= min, (2.9) 

where v_ = pu_ is the mass flow and w is a weighting function to be chosen. 
In this setup the variables are the potential <{> and the mass flow v. 

The density in (2.9) 

p = p(|grad <f)|) 

is obtained from Bernoulli's equation, i.e., 

P Y_1 = |i - M 2 ( | grad 4> 1 2 " l)J- 

Thus, (2.9) is a nonlinear least squares formulation, which is appropriate 
since it reflects the nonlinear character of transonic flow. Once a grid is 
selected (specific examples are given in the next section), the minimization 
of (2.9) over the associated finite element space leads to a nonlinear system 


/ 

n 


( 

V 

2 ?) 


— - grad 4 > 

+ w div v } 

( 

P 

" ) 


K($)«> = F 


( 2 . 10 ) 



In all of the numerical examples reported in the next section, (2.10) was 
solved by a combination of Newton's method and elimination. Issues related to 
this choice for the equation solver will be discussed in the next section. 

There are three main cases that are considered in this paper: 

Case 1 : smooth subsonic flows, 

Case 2 : smooth transonic flows, 

Case 3 : transonic flows with shocks. 

In the first case (2.9) can be used without modification, and in 
particular no weighting function is needed (i.e., w = 1 can be used). One 
does need special grids to obtain optimal accuracy (see [1]), and the criss- 
cross grid pattern which satisfies the grid decomposition property of [1] is 
used. 

In the second case a hyperbolic region appears but the flow field remains 
smooth. In this case there is a loss of accuracy in the hyperbolic region. 
In particular, with linear elements the pointwise accuracy in the mass flow v_ 
drops from 0(h) — in a generic mesh spacing — to 0(h). This can be corrected 
with a suitable choice of weighting function w, and details are given in [8]. 
This modification was not used in the results reported in this paper since the 
hyperbolic regions in question were too small for the suboptimal accuracy to 
have a major effect on the qualitative features of the flow. 

The third case is, by a wide margin, the most important as well as the 
most challenging. Here we use a weight w so that the term 



f , l2 

u 

2 ) 

/ 1 
n | 

|w|div u[ + 

- - grad <|> 

1 


( 2 . 11 ) 
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remains meaningful. In addition, modification to the density p = p(|grad 4> | ) 
must be introduced so that nonphysical expansion shocks are eliminated. 

For the choice of the weight w, we follow the developments introduced in 
[7]. For most flows, v = p_u is continuous across the shock [10]. Neverthe- 
less, It does not follow that div v_ is square integrable, and the primary 
rule derived from [7] is that w be chosen so that 

/ w|div vj 2 < «. (2.12) 
0 

This requires that w vanishes appropriately on the shock, which in turn 
means that (2.11) is a least squares principle in a degenerate L norm. A 
point of significance, on the other hand, is the fact that if w vanishes to 

O 

minimal order on the shock (in that (2.12) still holds), then optimal 0(h) 
can be achieved in unweighted L norms provided appropriate mesh refinement 
is introduced. This has been proved rigorously only in special cases (see 
[7]), yet the numerical results in the next section seem to indicate that the 
principle is general. 

These modifications alone do not yield an accurate simulation of the flow 
problem. To do this one must deal with the presence of nonphysical expansion 
shocks. In effect, (2.9) does not have a unique minimum, neither over 
infinite-dimensional function spaces nor over the finite-dimensional finite- 
element spaces. One can have expansion shocks, compression shocks, or both. 
What is interesting is the results In the next section tend to indicate that 
the case where both type of shocks appear tends to be the stable mode for 
(2.10). That is, an arbitrary choice of starting vector for Newton's methods 
applied to (2.10) tends to converge to this solution. 
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To eliminate expansion shocks we consider density biasing which in effect 
introduces streamwise diffusion into (2.1) - (2.2). Following [11] (see also 
[12] - [14]) the modified density takes the form 

p = p - pp g As, (2.13) 

where p g is the derivative of the density p along the streamwise 
direction. Since the density has the form 

p = p( | grad <j>|), 

the derivative p formally involves second derivatives of <f» • Since <j> is 
s 

expanded in terms of linear elements, it is necessary to replace p with a 
streamwise difference quotient; i.e., 


p = p - y Ap As, 


in the least squares formulation. 


(2.13') 


3. NUMERICAL RESULTS 

To illustrate the above ideas we selected the classic problem of a planar 
flow past a cylinder. The flow region plus boundary conditions are given in 
Figure 3.1. The configuration shown in this figure assumes that both the 
outflow and inflow remain subsonic. Figure 3.2 contains a typical grid. For 
economy only the top part of the flow region is shown, and the special 


refinement needed for the shocks is not shown. 
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The first set of results shows a typical subsonic flow pattern. The 

results are given in Figure 3.3 for a free stream Mach number of 

M =0.1. 

00 

Convergence studies at such Mach numbers are reported in [5] - [6]. These 

results indicate, with the type of grid shown in Figure 3.2, one can readily 

o 

achieve L error of 1% or less for the velocity field. 

The next set of results deal with the smooth transonic case. Of special 
interest here is the ability of the scheme to detect the onset of supersonic 
flow. Analytical techniques (see [15] and [16]) have given accurate values 

for the critical free stream Mach number M* as a function of d/D, where 

d is the diameter of the cylinder and D is the width of the channel. These 
results are reproduced in Figure 3.4. Numerical results from the least 
squares scheme are given in Figures 3.5 - 3.7 for M^ = .42, .45, and .50, 
respectively. The d/D ratio used for this case is 1/6. Extrapolation 
based on these results indicates that the critical Mach number is 
approximately .41, which is good agreement with Figure 3.4. 

The next set of results show what least squares based schemes produce when 
diffusion via density modification is not used. These are shown in Figure 3.8 
which contains plots of the velocity q = |u| versus angle 6 along the 
cylinder and at a radius slightly above the cylinder. The free stream Mach 
number is M^ = .5. The shock at the front of the cylinder is an expansion 
shock and is nonphysical. The one at the rear is a compression shock. A 
remarkable feature of this approximation is that the physically relevant 
compression shock is approximately in its correct position and is apparently 
unaffected by, the spurious shock. (Compare Figures 3.8 and 3.9.) 
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The solution shown in Figure 3.8 is apparently a stable mode for the 
nonlinear system (2.10). Indeed, Newton's method converged to this solution 
rather rapidly for a wide variety of initial conditions. 

In this regard, it is interesting to note that for the least squares 
formulation the Jacobian is not singular near the solution shown in Figure 
3.8. Density modifications are needed to remove the spurious shock shown at 
the front of the cylinder. However, they are not needed to obtain nonsingular 
Jacobians. 

The final results deal with the complete least squares system with the 
density modification discussed in the previous section. Figures 3.9 - 3.11 
show the velocity field over the cylinder, at a radius slightly larger that 
that of the cylinder, and at a radius in the free stream. Note that the 
spurious expansion shock has been totally eliminated. Moreover, the shock 
location and strength as well as the velocity profile appear to be correct as 
is the supersonic bubble shown in Figure 3.12. 

While we regard these numerical experiments as successful, there are a 
number of areas where the approach could be improved. The first issue 
concerns the equation solver. Once the density modification were introduced, 
the number of iterations increased by a factor of 2 to 3. Moreover, the 
solution shown in Figure 3.9 tended to be less "attractive" to the Newton 
iterations than that shown in Figure 3.8 (without density modifications). In 
fact, it was not difficult to find starting vectors where nonconvergence was 
seen, in the former case, although the starting state of a uniform flow always 
leads to convergence. This suggests that an alternative equation solver 
(e.g., preconditioned conjugate gradient) might be a more efficient choice for 
the equation solver. 
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A second issue concerns post-shock oscillations. These are seen in Figure 
3.10, which is the radius where the oscillations were found to be the most 
significant. These oscillations were not seen on the body of the cylinder 
(Figure 3.9) and disappeared rather rapidly away from the cylinder (Figure 
3.11). This is clearly a grid effect due to the slight misalignment of shock 
and grid. 


4. CONCLUSIONS 

Finite difference approximations to transonic flow problems are well- 
developed and have been successfully used for a wide range of problems. 

Nevertheless, there is still a need to develop finite element approaches for 
such problems for a variety of applications. We feel that the results 
presented here do show that such schemes can give physically meaningful 

simulations. 

On the other hand, our experience has tended to indicate that 

straightforward application of the basic finite element idea may not always be 
successful. Key computational issues are as follows: 

(i) There is a need to carefully develop the spaces in which the 

2 

approximations are formulated. Classical L spaces are generally 
inappropriate. 

(ii) Some form of diffusion (via density modifications or otherwise) appears 
to be needed. Moreover, care is needed in the way this diffusion is 
introduced. 

(iii) The geometrical pattern of the grid selected is of importance. Some 



- 11 - 


patterns are definitely superior to others. 

Finally, there are some important "bottlenecks" associated with the scheme 
employed in this paper, which, if properly addressed, could lead to an even 
more efficient approach. These include the following: 

(i) There is a need for an equation solver that is more efficient than the 
Newton method used in this paper. 

(ii) There is a need for adaptive grid refinement techniques that would lead 
to a better shock grid alignment than that achieved in this paper. 
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Figure 3.2. 


512 elements, 281 nodes, 


h = 0.30907 x 10 1 . 
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Figure 3.5. Plots of the supersonic pocket for = 0.42. 
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Figure 3.8. Velocity as a function of angle: (a) on cylinder, (b) 

slightly off cylinder — M = .51. 
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Figure 3.9. 
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squares scheme with density modification — = .5. 
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Figure 3.10. Velocity as a function of angle slightly off cylinder — full 
least squares scheme with density modification — M = .5. 
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Figure 3.11. Velocity as a function of angle half radius above cylinder 
full least squares scheme with density modification — = 
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